###########################################################################################
# CJR_on_UN.R
# Author: G Rosas
# January 2013
# Purpose: 
# Runs the CJR IRT model, along with different versions of competing principals,
# on data from votes in the UN General Assembly in 1970
# Though every step can be recreated running this file, the longer JAGS runs have been
# saved in a handful of RData files:
# UN2.RData contains the roll calls prepared by Y Shomer
# UNfixBetaRun.RData contains results from the CJR IRT model
# UNCompeteRun.RData is a basic competing principals model
# UNCompeteRunDiff.RData adds principal-specific probs of voting given agree/disagree
# UNCompeteRunDiffRecursive.RData adds self-actualizing vote choices for principals
#         (only for missing values)
###########################################################################################

# The original data for this section is in alliance.csv and comes from:
# Curtis Signorino and Jeffery Ritter. 1999.
# "Tau-b or Not Tau-b: Measuring the Similarity of Foreign Policy Positions",
# International Studies Quarterly, 43, pp. 115-144. 

# The data has been re-codified and it appears in "UN alliance data SUG.csv"
# To re-codify these data in a suitable format, we proceeded as follows:

# From the dyadic data in alliance.csv, we create four variables
# that code a country's status vis-a-vis the US and the USSR. 

# 1: defense alliance 
# 2: neutrality pact
# 3: alliance with no intervention
# 4: no alliance


###########################
# CJR on UN for 1970
###########################

# Require packages
library (MCMCpack)
library (MASS)
library (FastImputation)
library (superdiag)
library (coda)
library (mcmcplots)
library (sm)
library (foreign)
library (rjags)
library (runjags)
library (multicore)
library (random)
library (car)

load("UN2.RData")

rc <- Data.1970[,grep("V11", colnames(Data.1970))]
colnames(rc) <- paste("V", 1:ncol(rc), sep="")
rownames(rc) <- c(UnfactorColumns(Data.1970)$NATION[-127], "FIJI")
rc <- rc[-grep("MALDIVE.ISLANDS", rownames(rc)),]

codes    <- read.table ("UN_codebook.txt", header=TRUE, sep="\t")
alliance <- read.csv ("UN alliance data SUG.csv", header=TRUE)

# Build an expanded roll call matrix with country and alliance codes
rc.exp <- data.frame (rc, country=rownames(rc)); rc.exp$country <- as.character (rc.exp$country)
rownames (rc.exp) <- 1:nrow(rc.exp)

# Check names
cbind (rc.exp$country, rownames (rc))

# Do we have codes for all countries?
is.element (rc.exp$country, codes$NATION)
rc.exp <- merge (rc.exp, codes, by.x="country", by.y="NATION", all.x=TRUE, sort=T) # Notice that result is sorted
rc.exp <- merge (rc.exp, alliance, by.x="CODE", by.y="code", all.x=TRUE, sort=F)
rc.exp$US.client <- ifelse (!is.na(rc.exp$AllianceUS) & rc.exp$AllianceUS==1, 1, 0) 
rc.exp$Soviet.client <- ifelse (!is.na(rc.exp$AllianceSoviet) & rc.exp$AllianceSoviet==1, 1, 0) 
# Check names again

# Aye votes/percentage
Aye.pc <- rowSums (rc, na.rm=T) / rowSums (!is.na(rc))
# Abstention percentage
Abs.pc <- rowSums (is.na(rc), na.rm=T) / ncol (rc)

# Uncomment to plot abstention rates of countries
# plot (Abs.pc~Aye.pc, type="n", xlab="% Aye", ylab="% Abstention")
# text (xy.coords (Aye.pc, Abs.pc), labels=tolower(rownames (rc)))

################################
# Find potential anchor votes
################################
# We look for votes that pit one majority of US.party against a majority of Soviet.party
# These votes will anchor the spatial model

US.party <- length (rc.exp$US.client[rc.exp$US.client==1])
Soviet.party <- length (rc.exp$Soviet.client[rc.exp$Soviet.client==1])

mean.US.party.rate  <- colSums(rc.exp[rc.exp$US.client==1,c(3:70)], na.rm=T)/US.party
mean.Sov.party.rate <- colSums(rc.exp[rc.exp$Soviet.client==1,c(3:70)], na.rm=T)/Soviet.party
cbind (mean.US.party.rate, mean.Sov.party.rate)[c(7,11,13,22,23,24,57)]  
rc.exp$V1[rc.exp$Soviet.client==1] # Vote 1 is a firm candidate
rc.exp$V1[rc.exp$US.client==1] #
rc.exp$V1[rc.exp$country=="U.S.S.R."]
rc.exp$V1[rc.exp$country=="UNITED.STATES"]

rc.exp$V57[rc.exp$Soviet.client==1] # Vote 57 is the second firm candidate
rc.exp$V57[rc.exp$US.client==1] #
rc.exp$V57[rc.exp$country=="U.S.S.R."]
rc.exp$V57[rc.exp$country=="UNITED.STATES"]

# The following line is used much later, around line 700, for Voeten model
y.prelim <- rc.exp[,3:70]
Voeten.country.order <- rc.exp$country


######################################
# Run CJR model using MCMCpack
######################################

# This is a dry run using MCMCpack to use as a baseline
# The model can be re-run by commenting out the following lines of code
# We also include (commented out) code to assess convergence
# and to plot ideal points of countries

# result.1 <- MCMCirt1d(rc,
#   theta.constraints=list(U.S.S.R.=-2, UNITED.STATES=2),
#   verbose=10000,
#   T0=.25,
#   thin=50,
#   burnin=50000,
#   mcmc=100000, seed=197)
# 
# result.2 <- MCMCirt1d(rc,
#   theta.constraints=list(U.S.S.R.=-2, UNITED.STATES=2),
#   verbose=10000,
#   T0=.25,
#   thin=50,
#   burnin=50000,
#   mcmc=100000, seed=1971)

## THIS IS WHERE ALL OF THE DIFFERENT CHAINS SHOULD BE LOADED
# chainsConv <- mcmc.list(list (result.1[,-c(1,47)], result.2[,-c(1,47)]))  
# MixChains <- rbind (chainsConv[[1]], chainsConv[[2]])

# gelman.diag (chainsConv) # Good convergence, low values of Gelman-Rubin hat statistic for each parameter.  Multivariate psrf at around 1.3
# geweke.diag (MixChains)
# gw <- geweke.diag (MixChains)[[1]]
# mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 0.04

## Plot of Geweke stats
# hist (gw, 30, prob=TRUE)
# lines (density(gw))
# qqnorm (gw)

## Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUN"   # <= Change name here
#           , random=10)

# Theta.sig.irt <- matrix (NA, ncol=3, nrow=ncol(MixChains))
# for (i in 1:ncol(MixChains)){
#   Theta.sig.irt[i,] <- quantile (MixChains[,i], prob=c(0.25,0.5,0.75))
# }

# ord.ideal <- order (Theta.sig.irt[,2])
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,124), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.sig.irt[ord.ideal,1], x1=Theta.sig.irt[ord.ideal,3], y0=1:124, y1=1:124)
# points ( xy.coords ( Theta.sig.irt[ord.ideal,2], 1:124), pch=19)
# segments (x0=Theta.sig.irt[ord.ideal,1], x1=Theta.sig.irt[ord.ideal,3], y0=1:124, y1=1:124)
# text (  xy.coords ( Theta.sig.irt[ord.ideal,2]+1, 1:124), label=tolower(rownames(rc)[-c(1,47)])[ord.ideal], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")


######################################################################
###  Modeling abstentions with complete-data comp. principals model
######################################################################

# The previous dry run with MCMCpack looks good, with very extreme positions for very close allies
# We now need to see if we can improve on what we have
# One problem is that we need to fix scale based on ITEMS, not COUNTRIES, as we did before

# First, we will re-arrange voting data so that Soviet Clients appear last
rc.exp <- rc.exp[with(rc.exp, order(Soviet.client)), ]

# This is the regular IRT model, run again for comparison, using exact same priors
irt="model { 
  for (i in 1:n.legislators) {
    for (j in 1:n.item) {
      rc[i,j] ~ dbern(p[i,j]);
      probit(p[i,j]) <- beta[j]*theta[i] - alpha[j];
    }
  }
# PRIORS
  for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
  beta[fix.bill[1]] <- 4
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
  beta[fix.bill[2]] <- -4
  for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
  for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"

n.legislators <- nrow (rc.exp[,3:70])
n.item <- ncol (rc.exp[,3:70])
fix.bill <- c(1,57)

IRT.country.order <- as.character (rc.exp$country)

irt.parameters = c("theta", "alpha", "beta", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
irt.data <- dump.format(list(rc=as.matrix(rc.exp[,3:70]), n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill))

## This next step fits the model, and it is time-consuming. Instead of re-running,
## simply import the RData file named "UNfixBetaRun.RData", which has object "out"

# system.time(
#   out <- mclapply(1:2, function(x) {
#     jags.irt <- try(run.jags(   
#       model=irt,
#       monitor=irt.parameters, n.chains=1,
#       data=irt.data,
#       inits=irt.inits(),
#       thin=100, burnin=50000, sample=100000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.irt,"try-error")) {return()}
#     return(jags.irt)
#   }, mc.cores=2 )
# )
# save (out, file="UNfixBetaRun.RData")
load("UNfixBetaRun.RData")

chainsConv <- mcmc.list(list (out[[1]][[1]][[1]][,-c(fix.bill+194)], out[[2]][[1]][[1]][,-c(fix.bill+194)]))  
MixChains <- rbind (out[[1]][[1]][[1]], out[[2]][[1]][[1]])

## Convergence diagnostics
gelman.diag (chainsConv, multivariate=FALSE)
gw <- geweke.diag (MixChains)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 

# Plot of Geweke stats
# hist (gw, 30, prob=TRUE)
# lines (density(gw, na.rm=T))
# qqnorm (gw)

## Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUNfixBeta"   # <= Change name here
#           , random=10)
# Could use more thinning, but chains are mixing well


Theta <- MixChains[,c(1:126)]
Theta.irt <- matrix (NA, ncol=3, nrow=ncol(Theta))
for (i in 1:ncol(Theta)){
  Theta.irt[i,] <- quantile (Theta[,i], prob=c(0.25,0.5,0.75))
}

ord.ideal <- order (Theta.irt[,2])
# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.irt[ord.ideal,1], x1=Theta.irt[ord.ideal,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.irt[ord.ideal,2], 1:126), pch=19)
# segments (x0=Theta.irt[ord.ideal,1], x1=Theta.irt[ord.ideal,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.irt[ord.ideal,2]+0.6, 1:126), label=tolower(IRT.country.order[ord.ideal]), cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")


###########################################################################
# Competing principals model
###########################################################################

# This is the model with beta-distributed priors on the deltas
compete="model {
  for (i in 1:n.legislators) {
    for (j in 1:n.item) {
      probit(CJR[i,j]) <- beta[j]*theta[i] - alpha[j];
    }
  }
  for (k in 1:2) {
    for (j in 1:n.item) {
      lead[k,j] ~ dbern (0.5)
    }
  }
  for (i in 1:n.partisans) {
    for (j in 1:n.item) {
      Y[i,j] ~ dcat(Q[i,j,1:3]);         # code Nay as 1, NA as 2, Aye as 3
      Q[i,j,1] <- pow(probVoteDisagree*(1-CJR[i,j]), p[i,j]) * pow(probVoteAgree*(1-CJR[i,j]), (1 - p[i,j]));
      Q[i,j,2] <- 1 - Q[i,j,1] - Q[i,j,3];
      Q[i,j,3] <- pow(probVoteAgree*CJR[i,j], p[i,j]) * pow(probVoteDisagree*CJR[i,j], (1 - p[i,j]));
      p[i,j] <- lead[party.membership[i],j]; 
    }
  }
  for (i in (n.partisans+1):n.legislators) {
    for (j in 1:n.item) {
      Y[i,j] ~ dbern(CJR[i,j]);
    }
  }
# priors
#   delta[1] ~ dbeta(1,5); # This prior puts a lot of mass on values close to 0
#   delta[2] ~ dbeta(5,1); # This prior puts a lot of mass on values close to 1
delta[1] ~ dbeta(1,1); 
delta[2] ~ dbeta(1,1); 
probVoteDisagree <- delta[1];
probVoteAgree    <- delta[2];
for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
  beta[fix.bill[1]] <-  4
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
  beta[fix.bill[2]] <- -4
  for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"

party.membership <- ifelse (rc.exp$US.client==1, 1, ifelse(rc.exp$Soviet.client==1, 2, 3)) 
new.ord <- order (party.membership)
party.membership <- party.membership[new.ord]
RC <- rc.exp[new.ord,3:70]
rownames (RC) <- rc.exp$country[new.ord]
new.RC <- RC #  RC[-c(40,50),]  # Use this if USSR (50) and US (40) are omitted

#Quick checks
party.membership[rownames(new.RC)=="UNITED.STATES"]  # should be 1
party.membership[rownames(new.RC)=="U.S.S.R."]       # should be 2

lead.SOV <- RC["U.S.S.R.",]
lead.USA <- RC["UNITED.STATES",]
lead <- rbind(lead.USA, lead.SOV)
# code NA as 2, Nay as 1, Aye as 3 (40 and 50 are rows for US and USSR in the re-arranged roll-call matrix)
y <- ifelse(is.na(new.RC)==TRUE, 2, ifelse(new.RC==0,1,3))    
n.legislators <- nrow (y)

# Use the next line if want to omit US and USSR as "whips"
# party.membership.no.whips <- party.membership[-c(19,52)]
party.membership.no.whips <- party.membership
n.partisans <- length (party.membership.no.whips[party.membership.no.whips<3])
n.item <- ncol (y)
Y1 <- y[1:n.partisans,]
Y2 <- new.RC[(n.partisans+1):n.legislators,]
Y  <- rbind (Y1, Y2)
fix.bill <- c(1,57)

compete.parameters = c("theta", "alpha", "beta", "probVoteAgree", "lead", "probVoteDisagree", "deviance")

compete.inits <- function() {
  dump.format(
    list(
        theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , delta = sort(runif (2)) 
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}

compete.data <- dump.format(list(Y=as.matrix(Y), lead=as.matrix(lead), party.membership=party.membership.no.whips[1:n.partisans], n.partisans=n.partisans, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill))

#  This next step already done, simply import the RData file named "UNCompeteRun.RData", which has object "out2"
# system.time(
#   out2 <- mclapply(1:2, function(x) {
#     jags.compete <- try(run.jags(   
#       model=compete,
#       monitor=compete.parameters,
#       n.chains=2,
#       data=compete.data,
#       inits=compete.inits(),
#       thin=100, burnin=50000, sample=100000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.compete,"try-error")) {return()}
#     return(jags.compete)
#   }, mc.cores=2 )
# )
# 
# save (out2, file="UNCompeteRun.RData")
load ("UNCompeteRun.RData")

# Do not include "lead", as there is no variation here
chain1 <- out2[[1]][[1]][[1]][,-grep("lead", colnames(out2[[1]][[1]][[1]]))]
chain2 <- out2[[2]][[1]][[1]][,-grep("lead", colnames(out2[[2]][[1]][[1]]))]
chainsConv2 <- mcmc.list(list (chain1[,-c(fix.bill+194)], chain2[,-c(fix.bill+194)]))  
MixChains2  <- rbind (chain1, chain2)

gelman.diag (chainsConv2) # This is still not working, figure out what's wrong
geweke.diag (MixChains2)
gw <- geweke.diag (MixChains2)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 0.057

# Plot of Geweke stats
# hist (gw, 30, prob=TRUE)
# lines (density(gw, na.rm=T))
# qqnorm (gw)

# Traceplots of theta parameters
# mcmcplot (chainsConv2, dir="~/Dropbox/UN/"   
# 		  , filename="traceplotsUNcompete"   # <= Change name here
# 		  , random=10)
# # Could use more thinning, but they are mixing well


Theta <- MixChains2[,c(1:126)]
Theta.compete <- matrix (NA, ncol=3, nrow=ncol(Theta))
for (i in 1:ncol(Theta)){
	Theta.compete[i,] <- quantile (Theta[,i], prob=c(0.25,0.5,0.75))
}

ord.compete <- order (Theta.compete[,2])

# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.compete[ord.compete,1], x1=Theta.compete[ord.compete,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.compete[ord.compete,2], 1:126), pch=19)
# segments (x0=Theta.compete[ord.compete,1], x1=Theta.compete[ord.compete,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.compete[ord.compete,2]+0.6, 1:126), label=tolower(rownames(new.RC))[ord.compete], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")
# 
# cor (cbind (Theta.compete[,2], Theta.irt[new.ord,2]))
# par (mar=c(3,3,1,1))
# plot (c(-4,4),c(-4,4),type="n")
# segments (x0=Theta.compete[,1], x1=Theta.compete[,3], y0=Theta.irt[new.ord,2], y1=Theta.irt[new.ord,2])
# segments (x0=Theta.compete[,2], x1=Theta.compete[,2], y0=Theta.irt[new.ord,1], y1=Theta.irt[new.ord,3])
# text (xy.coords(Theta.compete[,2], Theta.irt[new.ord,2]), label=tolower(rownames(new.RC)), cex=0.6)

###########################################################################
# Competing principals model, with whip-specific probabilities of voting
###########################################################################

# This is the model with beta-distributed priors on the deltas and different deltas per whip
compete.bis="model {
for (i in 1:n.legislators) {
for (j in 1:n.item) {
probit(CJR[i,j]) <- beta[j]*theta[i] - alpha[j];
}
}
for (k in 1:2) {
for (j in 1:n.item) {
lead[k,j] ~ dbern (0.5)
}
}
for (i in 1:n.partisans) {
for (j in 1:n.item) {
Y[i,j] ~ dcat(Q[i,j,1:3]);         # code Nay as 1, NA as 2, Aye as 3
Q[i,j,1] <- pow(probVoteDisagree[party.membership[i]]*(1-CJR[i,j]), p[i,j]) * pow(probVoteAgree[party.membership[i]]*(1-CJR[i,j]), (1 - p[i,j]));
Q[i,j,2] <- 1 - Q[i,j,1] - Q[i,j,3];
Q[i,j,3] <- pow(probVoteAgree[party.membership[i]]*CJR[i,j], p[i,j]) * pow(probVoteDisagree[party.membership[i]]*CJR[i,j], (1 - p[i,j]));
p[i,j] <- lead[party.membership[i],j]; 
}
}
for (i in (n.partisans+1):n.legislators) {
for (j in 1:n.item) {
Y[i,j] ~ dbern(CJR[i,j]);
}
}
# priors
#   delta[1] ~ dbeta(1,5); # This prior puts a lot of mass on values close to 0
#   delta[2] ~ dbeta(5,1); # This prior puts a lot of mass on values close to 1
  delta[1] ~ dbeta(1,1); 
  delta[2] ~ dbeta(1,1); 
  delta[3] ~ dbeta(1,1); 
  delta[4] ~ dbeta(1,1); 
  probVoteDisagree[1] <- delta[1];
  probVoteAgree[1]    <- delta[2];
  probVoteDisagree[2] <- delta[3];
  probVoteAgree[2]    <- delta[4];
for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
beta[fix.bill[1]] <-  4
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
beta[fix.bill[2]] <- -4
for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"


compete.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , delta = sort(runif (4)) 
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}

## This next step already done, simply import the RData file named "UNCompeteRunDiff.RData", which has object "out2"
# system.time(
#   out3 <- mclapply(1:2, function(x) {
#     jags.compete <- try(run.jags(   
#       model=compete.bis,
#       monitor=compete.parameters,
#       n.chains=1,
#       data=compete.data,
#       inits=compete.inits(),
#       thin=100, burnin=50000, sample=100000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.compete,"try-error")) {return()}
#     return(jags.compete)
#   }, mc.cores=2 )
# )
# 
# save (out3, file="UNCompeteRunDiff.RData")

load ("UNCompeteRunDiff.RData")

chain1 <- out3[[1]][[1]][[1]][,-grep("lead", colnames(out3[[1]][[1]][[1]]))]
chain2 <- out3[[2]][[1]][[1]][,-grep("lead", colnames(out3[[2]][[1]][[1]]))]

Lead1.partial <-  out3[[1]][[1]][[1]][,grep("lead", colnames(out3[[1]][[1]][[1]]))]
Lead2.partial <-  out3[[2]][[1]][[1]][,grep("lead", colnames(out3[[2]][[1]][[1]]))]

chainsConv3 <- mcmc.list(list (chain1[,-c(fix.bill+194)], chain2[,-c(fix.bill+194)]))  
MixChains3  <- rbind (chain1, chain2)

gelman.diag (chainsConv3)
gw <- geweke.diag (MixChains3)[[1]]
mean(abs(gw)>2, na.rm=TRUE) 
# should be close to 0.05; it is 0.071, so could improve

# Plot of Geweke stats
# hist (gw, 30, prob=TRUE)
# lines (density(gw, na.rm=T))
# qqnorm (gw)

# Traceplots of theta parameters
# mcmcplot (chainsConv3, dir="~/Dropbox/UN/"   
# 		  , filename="traceplotsUNcompeteDiff"   # <= Change name here
# 		  , random=10)
# # Could use more thinning, but they are mixing well


Theta <- MixChains3[,c(1:126)]
Theta.compete.bis <- matrix (NA, ncol=3, nrow=ncol(Theta))
for (i in 1:ncol(Theta)){
	Theta.compete.bis[i,] <- quantile (Theta[,i], prob=c(0.25,0.5,0.75))
}

ord.compete.bis <- order (Theta.compete.bis[,2])
# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.compete.bis[ord.compete.bis,1], x1=Theta.compete.bis[ord.compete.bis,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.compete.bis[ord.compete.bis,2], 1:126), pch=19)
# segments (x0=Theta.compete.bis[ord.compete.bis,1], x1=Theta.compete.bis[ord.compete.bis,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.compete.bis[ord.compete.bis,2]+0.6, 1:126), label=tolower(rownames(new.RC))[ord.compete.bis], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")
# 
# 
# cor (cbind (Theta.compete.bis[,2], Theta.compete[,2], Theta.irt[new.ord,2]))
# par (mar=c(3,3,1,1))
# plot (c(-4,4),c(-4,4),type="n")
# segments (x0=Theta.compete.bis[,1], x1=Theta.compete.bis[,3], y0=Theta.irt[new.ord,2], y1=Theta.irt[new.ord,2])
# segments (x0=Theta.compete.bis[,2], x1=Theta.compete.bis[,2], y0=Theta.irt[new.ord,1], y1=Theta.irt[new.ord,3])
# text (xy.coords(Theta.compete.bis[,2], Theta.irt[new.ord,2]), label=tolower(rownames(new.RC)), cex=0.6)
# 
# plot (c(-4,4),c(-4,4),type="n")
# segments (x0=Theta.compete[,1], x1=Theta.compete[,3], y0=Theta.compete.bis[,2], y1=Theta.compete.bis[,2])
# segments (x0=Theta.compete[,2], x1=Theta.compete[,2], y0=Theta.compete.bis[,1], y1=Theta.compete.bis[,3])
# text (xy.coords(Theta.compete[,2], Theta.compete.bis[,2]), label=tolower(rownames(new.RC)), cex=0.6)



###########################################################################
# Competing principals model, with whip-specific probabilities of voting
# and with "lead" values conditional on Y values
###########################################################################


# This is the model with beta priors on the deltas, different deltas per whip, and self-actualizing values for missing "lead"
compete.3="model {
  for (i in 1:n.legislators) {
  for (j in 1:n.item) {
    probit(CJR[i,j]) <- beta[j]*theta[i] - alpha[j];
    }
  }
  for (k in 1:2) {
    for (j in 1:n.item) {
      lead[k,j] ~ dbern (CJR[principal[k],j])
    }
  }
  for (i in 1:n.partisans) {
    for (j in 1:n.item) {
      Y[i,j] ~ dcat(Q[i,j,1:3]);         # code Nay as 1, NA as 2, Aye as 3
      Q[i,j,1] <- pow(probVoteDisagree[party.membership[i]]*(1-CJR[i,j]), p[i,j]) * pow(probVoteAgree[party.membership[i]]*(1-CJR[i,j]), (1 - p[i,j]));
      Q[i,j,2] <- 1 - Q[i,j,1] - Q[i,j,3];
      Q[i,j,3] <- pow(probVoteAgree[party.membership[i]]*CJR[i,j], p[i,j]) * pow(probVoteDisagree[party.membership[i]]*CJR[i,j], (1 - p[i,j]));
      p[i,j] <- lead[party.membership[i],j]; 
    }
  }
  for (i in (n.partisans+1):n.legislators) {
    for (j in 1:n.item) {
      Y[i,j] ~ dbern(CJR[i,j]);
    }
  }
# priors
  for (m in 1:4) {
    delta[m] ~ dbeta(1,1);
  }
  probVoteDisagree[1] <- delta[1];
  probVoteAgree[1]    <- delta[2];
  probVoteDisagree[2] <- delta[3];
  probVoteAgree[2]    <- delta[4];
  for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
  beta[fix.bill[1]] <-  4
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
  beta[fix.bill[2]] <- -4
  for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
  for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"


compete.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , delta = sort(runif (4)) 
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}

compete.data <- dump.format(list(Y=as.matrix(Y), lead=as.matrix(lead), principal=c(40,50), party.membership=party.membership.no.whips[1:n.partisans], n.partisans=n.partisans, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill))
compete.parameters = c("theta", "alpha", "beta", "probVoteAgree", "lead", "probVoteDisagree", "deviance")

#  This next step already done, simply import the RData file named "UNCompeteRunDiffRecursive.RData", which has object "out2"
# system.time(
#   out4 <- mclapply(1:2, function(x) {
#     jags.compete <- try(run.jags(   
#       model=compete.3,
#       monitor=compete.parameters,
#       n.chains=1,
#       data=compete.data,
#       inits=compete.inits(),
#       thin=100, burnin=50000, sample=100000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.compete,"try-error")) {return()}
#     return(jags.compete)
#   }, mc.cores=2 )
# )
# 
# save (out4, file="UNCompeteRunDiffRecursive.RData")
load ("UNCompeteRunDiffRecursive.RData")

chain1 <- out4[[1]][[1]][[1]][,-grep("lead", colnames(out4[[1]][[1]][[1]]))]
chain2 <- out4[[2]][[1]][[1]][,-grep("lead", colnames(out4[[2]][[1]][[1]]))]

Lead1.full <-  out4[[1]][[1]][[1]][,grep("lead", colnames(out4[[1]][[1]][[1]]))]
Lead2.full <-  out4[[2]][[1]][[1]][,grep("lead", colnames(out4[[2]][[1]][[1]]))]

chainsConv4 <- mcmc.list(list (chain1[,-c(fix.bill+194)], chain2[,-c(fix.bill+194)]))  
MixChains4  <- rbind (chain1, chain2)

gelman.diag (chainsConv4)
gw <- geweke.diag (MixChains4)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05, it is 0.038

# Uncomment to plot Geweke stats
# hist (gw, 30, prob=TRUE)
# lines (density(gw, na.rm=T))
# qqnorm (gw)

# Traceplots of theta parameters
# mcmcplot (chainsConv4, dir="~/Dropbox/UN/"   
#   	  , filename="traceplotsUNcompeteDiffRecursive"   # <= Change name here
# 		  , random=10)
# # Could use many more scans and lots more thinning, but they are mixing well; keep for the time being


Theta <- MixChains4[,c(1:126)]
Theta.compete.4 <- matrix (NA, ncol=3, nrow=ncol(Theta))
for (i in 1:ncol(Theta)){
  Theta.compete.4[i,] <- quantile (Theta[,i], prob=c(0.25,0.5,0.75))
}

ord.compete.4 <- order (Theta.compete.4[,2])
# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.compete.4[ord.compete.4,1], x1=Theta.compete.4[ord.compete.4,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.compete.4[ord.compete.4,2], 1:126), pch=19)
# text (  xy.coords ( Theta.compete.4[ord.compete.4,2]+0.6, 1:126), label=tolower(rownames(new.RC))[ord.compete.4], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")


# cor (cbind (Theta.compete.4[,2], Theta.compete.bis[,2], Theta.compete[,2], Theta.irt[new.ord,2]))
# par (mar=c(3,3,1,1))
# plot (c(-4,4),c(-4,4),type="n")
# segments (x0=Theta.compete.4[,1], x1=Theta.compete.4[,3], y0=Theta.irt[new.ord,2], y1=Theta.irt[new.ord,2])
# segments (x0=Theta.compete.4[,2], x1=Theta.compete.4[,2], y0=Theta.irt[new.ord,1], y1=Theta.irt[new.ord,3])
# text (xy.coords(Theta.compete.4[,2], Theta.irt[new.ord,2]), label=tolower(rownames(new.RC)), cex=0.6)
# 
# plot (c(-4,4),c(-4,4),type="n")
# segments (x0=Theta.compete.4[,1], x1=Theta.compete.4[,3], y0=Theta.compete.bis[,2], y1=Theta.compete.bis[,2])
# segments (x0=Theta.compete.4[,2], x1=Theta.compete.4[,2], y0=Theta.compete.bis[,1], y1=Theta.compete.bis[,3])
# text (xy.coords(Theta.compete.4[,2], Theta.compete.4[,2]), label=tolower(rownames(new.RC)), cex=0.6)

###########################################################################################
# Additions for AJPS RNR, January 10, 2013
# GR
# Per reviewer's request, a reestimation based on Voeten's model
###########################################################################################

# This is our interpretation of Voeten's model
voeten="model {
  for (i in 1:n.legislators) {
    for (j in 1:n.item) {
      y[i,j] ~ dcat(p[i,j,1:3]);
      p[i,j,1] <- 1 - Q[i,j,1];
      p[i,j,2] <- Q[i,j,1] - Q[i,j,2];
      p[i,j,3] <- Q[i,j,2];
      for (i.cut in 1:n.cut){
        logit(Q[i,j,i.cut]) <- beta[j]*theta[i] + alpha[j] - tau[j,i.cut];
      }
    }
  }
  # PRIORS
  for(j in 1:n.item) { 
    alpha[j] ~ dnorm(0, 0.25); 
    for(i.cut in 1:n.cut){
      tau.unsorted[j,i.cut] ~ dnorm(alpha[j],0.01)
    }
      tau[j,1:n.cut] <- sort(tau.unsorted[j,1:n.cut]-mean(tau.unsorted[j,1:n.cut]))
  }
  beta[fix.bill[1]] <- 4
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
  beta[fix.bill[2]] <- -4
  for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
  for(i in 1:n.legislators) { theta[i] ~ dnorm(0,1); }
}"


y <- matrix (NA, nrow=nrow(y.prelim), ncol=ncol(y.prelim))
for (i in 1:nrow(y.prelim)){
  y[i,] <- recode (y.prelim[i,], "NA=2; 0=1; 1=3")
}
n.legislators <- nrow (y)
n.item <- ncol (y)
fix.bill <- c(1,57)

irt.parameters = c("theta", "alpha", "beta", "tau", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , tau.unsorted = cbind (rnorm(n.item), rnorm(n.item))
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
irt.data <- dump.format(list(y=y, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill, n.cut=2))


#  This next step already done, simply import the RData file named "UNvoetenRun.RData", which has object "out2"
# system.time(
#   out <- mclapply(1:2, function(x) {
#     jags.irt <- try(run.jags(   
#       model=voeten,
#       monitor=irt.parameters, n.chains=1,
#       data=irt.data,
#       inits=irt.inits(),
#       thin=100, burnin=50000, sample=100000,
#       # thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.irt,"try-error")) {return()}
#     return(jags.irt)
#   }, mc.cores=2 )
# )
# save (out, file="UNvoetenRun.RData")
load("UNvoetenRun.RData")

chainsConv <- mcmc.list(list (out[[1]][[1]][[1]][,-c(fix.bill+194)], out[[2]][[1]][[1]][,-c(fix.bill+194)]))  
MixChains5 <- rbind (out[[1]][[1]][[1]], out[[2]][[1]][[1]])

gelman.diag (chainsConv, multivariate=FALSE) # should all be close to 1; they are
geweke.diag (MixChains5)
gw <- geweke.diag (MixChains5)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 

# # Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUNvoeten"   # <= Change name here
#           , random=10)
# Excellent mixing


Theta.voeten <- matrix (NA, ncol=3, nrow=ncol(MixChains5[,c(1:126)]))
for (i in 1:ncol(MixChains5[,c(1:126)])){
  Theta.voeten[i,] <- quantile (MixChains5[,c(1:126)][,i], prob=c(0.25,0.5,0.75))
}


ch.ord <- numeric()
for (i in 1:length(Voeten.country.order)) {
  cname <- IRT.country.order[i]
  ch.ord[i] <- which (Voeten.country.order==cname)
}

Voeten.country.order.2 <- Voeten.country.order[ch.ord]


ord.voeten <- order (Theta.voeten[,2])
# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# segments (x0=Theta.voeten[ord.voeten,1], x1=Theta.voeten[ord.voeten,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.voeten[ord.voeten,2]+0.6, 1:126), label=tolower(Voeten.country.order)[ord.voeten], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")
# 
# plot ( c(-6,6), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.irt[ord.ideal,1], x1=Theta.irt[ord.ideal,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.irt[ord.ideal,2], 1:126), pch=19, cex=0.5)
# text (  xy.coords ( Theta.irt[ord.ideal,2]+0.6, 1:126), label=tolower(IRT.country.order)[ord.ideal], cex=0.5 )
# segments (x0=Theta.voeten[ch.ord,1][ord.ideal], x1=Theta.voeten[ch.ord,3][ord.ideal], y0=1:126, y1=1:126, col="grey")
# points ( xy.coords ( Theta.voeten[ch.ord,2][ord.ideal], 1:126), pch=19, cex=0.5, col="gray")
# text (  xy.coords ( Theta.voeten[ch.ord,2][ord.ideal]-1.6, 1:126), label=tolower(Voeten.country.order.2)[ord.ideal], cex=0.5 )


####################################################################
# AJPS RNR
# Voeten model, second try
####################################################################
# Only abstentions coded exclusively as 2, all the rest are missing values
# Data_1970.1 is an object in "UN2.RData" that contains the original codes for GA voting in 1970
# We only need to convert 0, 5 and 6 into missing values.
# 2 remains as 2 (these are actual abstentions). Turn 1 into 3 (AYE) and 3 into 1 (NAY) to keep consistency
Data.Voeten <- as.matrix (Data_1970.1[,-1])
Data.Voeten <- recode (Data.Voeten, "1=3; 3=1; 2=2; else=NA")
# Data.Voeten$country.code <- Data_1970_final[1]  # recall order of countries, which should be identical to the one we used for regular IRT.

y <- Data.Voeten[-grep("MALDIVE.ISLANDS", c(UnfactorColumns(Data.1970)$NATION[-127], "FIJI")),]
y.country.order <- c(UnfactorColumns(Data.1970)$NATION[-127], "FIJI")[-grep("MALDIVE.ISLANDS", c(UnfactorColumns(Data.1970)$NATION[-127], "FIJI"))]
n.legislators <- nrow (y)
n.item <- ncol (y)
fix.bill <- c(1,57)

y.inits <- recode (y, "NA=2; else=NA")
irt.parameters = c("theta", "alpha", "beta", "tau", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , y = y.inits + matrix (rbinom(ncol(y.inits)*nrow(y.inits), 1, 0.5), ncol=ncol(y.inits), nrow=nrow(y.inits))
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , tau.unsorted = cbind (rnorm(n.item), rnorm(n.item))
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
irt.data <- dump.format(list(y=y, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill, n.cut=2))


#  This next step already done, simply import the RData file named "UNvoetenRunAbstentionsOnly.RData", which has object "out"
# system.time(
#   out <- mclapply(1:2, function(x) {
#     jags.irt <- try(run.jags(   
#       model=voeten,
#       monitor=irt.parameters, n.chains=1,
#       data=irt.data,
#       inits=irt.inits(), adapt=10000,
#       thin=50, burnin=50000, sample=1000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.irt,"try-error")) {return()}
#     return(jags.irt)
#   }, mc.cores=2 )
# )
# save (out, file="UNvoetenRunAbstentionsOnly.RData")

load ("UNvoetenRunAbstentionsOnly.RData")

chainsConv <- mcmc.list(list (out[[1]][[1]][[1]][,-c(fix.bill+194)], out[[2]][[1]][[1]][,-c(fix.bill+194)]))  
MixChains6 <- rbind (out[[1]][[1]][[1]], out[[2]][[1]][[1]])

gelman.diag (chainsConv, multivariate=FALSE) # should all be close to 1; they are
gw <- geweke.diag (MixChains6)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 0.0503

# Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUNvoeten2"   # <= Change name here
#           , random=10)
# Excellent mixing

Theta.voeten.2 <- matrix (NA, ncol=3, nrow=ncol(MixChains6[,c(1:126)]))
for (i in 1:ncol(MixChains6[,c(1:126)])){
  Theta.voeten.2[i,] <- quantile (MixChains6[,c(1:126)][,i], prob=c(0.25,0.5,0.75))
}

# Ordered preferences according to Voeten's model, based on abstentions only as 2s
ord.voeten.2 <- order (Theta.voeten.2[,2])

# To recover an order comparable to other runs
ch.ord.voeten.2 <- numeric()
for (i in 1:length(y.country.order)) {
  cname <- IRT.country.order[i]
  ch.ord.voeten.2[i] <- which (y.country.order==cname)
}

y.country.order.2 <- y.country.order[ch.ord.voeten.2]

# Some comparisons
# plot (Theta.voeten.2[ch.ord.voeten.2,2], Theta.voeten[ch.ord,2])
# plot (Theta.irt[,2], Theta.voeten[ch.ord,2])
# plot (Theta.irt[,2], Theta.voeten.2[ch.ord.voeten.2,2])

# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-9,9), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# segments (x0=Theta.voeten.2[ord.voeten.2,1], x1=Theta.voeten.2[ord.voeten.2,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.voeten.2[ord.voeten.2,2]+1, 1:126), label=tolower(y.country.order)[ord.voeten.2], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")
# # 
# plot ( c(-9,9), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.irt[ord.ideal,1], x1=Theta.irt[ord.ideal,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.irt[ord.ideal,2], 1:126), pch=19, cex=0.5)
# text (  xy.coords ( Theta.irt[ord.ideal,2]+0.6, 1:126), label=tolower(IRT.country.order)[ord.ideal], cex=0.5 )
# segments (x0=Theta.voeten.2[ch.ord.voeten.2,1][ord.ideal], x1=Theta.voeten.2[ch.ord.voeten.2,3][ord.ideal], y0=1:126, y1=1:126, col="grey")
# points ( xy.coords ( Theta.voeten.2[ch.ord.voeten.2,2][ord.ideal], 1:126), pch=19, cex=0.5, col="gray")
# text (  xy.coords ( Theta.voeten.2[ch.ord.voeten.2,2][ord.ideal]-1.6, 1:126), label=tolower(y.country.order.2)[ord.ideal], cex=0.5 )


# This model is the closest we can get to Voeten. It mixes well, but delivers results (theta positions) that are nowhere anything that makes sense
# We'll give this one last try based on tau intercepts that are fixed across votes (we will still accept a random alpha intercept)

voetenFixedTau="model {
  for (i in 1:n.legislators) {
    for (j in 1:n.item) {
y[i,j] ~ dcat(p[i,j,1:3]);
p[i,j,1] <- 1 - Q[i,j,1];
p[i,j,2] <- Q[i,j,1] - Q[i,j,2];
p[i,j,3] <- Q[i,j,2];
for (i.cut in 1:n.cut){
logit(Q[i,j,i.cut]) <- beta[j]*theta[i] + alpha[j] - tau[i.cut];
}
}
}
# PRIORS
for(j in 1:n.item) { 
  alpha[j] ~ dnorm(0, 0.25); 
}
for(i.cut in 1:n.cut){  
  tau.unsorted[i.cut] ~ dnorm(0,0.1)
}
tau[1:n.cut] <- sort(tau.unsorted[1:n.cut]-mean(tau.unsorted[1:n.cut]))
beta[fix.bill[1]] <- 4
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
beta[fix.bill[2]] <- -4
for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
for(i in 1:n.legislators) { theta[i] ~ dnorm(0,1); }
}"

y.inits <- recode (y, "NA=2; else=NA")
irt.parameters = c("theta", "alpha", "beta", "tau", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , y = y.inits + matrix (rbinom(ncol(y.inits)*nrow(y.inits), 1, 0.5), ncol=ncol(y.inits), nrow=nrow(y.inits))
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , tau.unsorted = rnorm(2)
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}


#  This next step already done, simply import the RData file named "UNvoetenRunAbstentionsOnlyFixedTau.RData", which has object "out5"
# system.time(
#   out5 <- mclapply(1:2, function(x) {
#     jags.irt <- try(run.jags(   
#       model=voetenFixedTau,
#       monitor=irt.parameters, n.chains=1,
#       data=irt.data,
#       inits=irt.inits(), adapt=10000,
#       thin=50, burnin=50000, sample=1000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.irt,"try-error")) {return()}
#     return(jags.irt)
#   }, mc.cores=2 )
# )
# save (out5, file="UNvoetenRunAbstentionsOnlyFixedTau.RData")

load ("UNvoetenRunAbstentionsOnlyFixedTau.RData")

chainsConv <- mcmc.list(list (out5[[1]][[1]][[1]][,-c(fix.bill+194)], out5[[2]][[1]][[1]][,-c(fix.bill+194)]))  
MixChains7 <- rbind (out5[[1]][[1]][[1]], out5[[2]][[1]][[1]])

gelman.diag (chainsConv, multivariate=FALSE) # should all be close to 1; they are
gw <- geweke.diag (MixChains7)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 0.0532

# Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUNvoetenFixedTau"   # <= Change name here
#           , random=10)
# Excellent mixing

Theta.voeten.fixedTau <- matrix (NA, ncol=3, nrow=ncol(MixChains7[,c(1:126)]))
for (i in 1:ncol(MixChains7[,c(1:126)])){
  Theta.voeten.fixedTau[i,] <- quantile (MixChains7[,c(1:126)][,i], prob=c(0.25,0.5,0.75))
}

# Ordered preferences according to Voeten's model, based on abstentions only as 2s
ord.voeten.tau <- order (Theta.voeten.fixedTau[,2])

# To recover an order comparable to other runs
ch.ord.voeten.tau <- numeric()
for (i in 1:length(y.country.order)) {
  cname <- IRT.country.order[i]
  ch.ord.voeten.tau[i] <- which (y.country.order==cname)
}

y.country.order.tau <- y.country.order[ch.ord.voeten.tau]

# Some comparisons
# plot (Theta.voeten.fixedTau[ch.ord.voeten.tau,2], Theta.voeten[ch.ord,2])
# plot (Theta.voeten.fixedTau[ch.ord.voeten.tau,2], Theta.voeten.2[ch.ord.voeten.2,2])
# plot (Theta.irt[,2], Theta.voeten.fixedTau[ch.ord.voeten.tau,2])

# Uncomment to plot
# par (mar=c(2,1,1,1), cex.axis=1, las=1)
# plot ( c(-9,9), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# segments (x0=Theta.voeten.fixedTau[ord.voeten.tau,1], x1=Theta.voeten.fixedTau[ord.voeten.tau,3], y0=1:126, y1=1:126)
# text (  xy.coords ( Theta.voeten.fixedTau[ord.voeten.tau,2]+1, 1:126), label=tolower(y.country.order)[ord.voeten.tau], cex=0.5 )
# abline (v=c(2,-2), lty=3, col="red")
# # 
# plot ( c(-9,9), c(1,126), type="n", ylab="", xlab="", main="", axes=F)
# axis(1)
# segments (x0=Theta.irt[ord.ideal,1], x1=Theta.irt[ord.ideal,3], y0=1:126, y1=1:126)
# points ( xy.coords ( Theta.irt[ord.ideal,2], 1:126), pch=19, cex=0.5)
# text (  xy.coords ( Theta.irt[ord.ideal,2]+0.6, 1:126), label=tolower(IRT.country.order)[ord.ideal], cex=0.5 )
# segments (x0=Theta.voeten.fixedTau[ch.ord.voeten.tau,1][ord.ideal], x1=Theta.voeten.fixedTau[ch.ord.voeten.tau,3][ord.ideal], y0=1:126, y1=1:126, col="grey")
# points ( xy.coords ( Theta.voeten.fixedTau[ch.ord.voeten.tau,2][ord.ideal], 1:126), pch=19, cex=0.5, col="gray")
# text (  xy.coords ( Theta.voeten.fixedTau[ch.ord.voeten.tau,2][ord.ideal]-1.6, 1:126), label=tolower(y.country.order.2)[ord.ideal], cex=0.5 )

# This last step returns a non-sensical order similar to the one in the previous exercise

y <- matrix (NA, nrow=nrow(y.prelim), ncol=ncol(y.prelim))
for (i in 1:nrow(y.prelim)){
  y[i,] <- recode (y.prelim[i,], "NA=2; 0=1; 1=3")
}
n.legislators <- nrow (y)
n.item <- ncol (y)
fix.bill <- c(1,57)

irt.parameters = c("theta", "alpha", "beta", "tau", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      , tau.unsorted = rnorm(2)
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
irt.data <- dump.format(list(y=y, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill, n.cut=2))


# This next step already done, simply import the RData file named "UNvoetenRunLastTry.RData", which has object "out5"
# system.time(
#   out6 <- mclapply(1:2, function(x) {
#     jags.irt <- try(run.jags(   
#       model=voetenFixedTau,
#       monitor=irt.parameters, n.chains=1,
#       data=irt.data,
#       inits=irt.inits(), adapt=10000,
#       thin=100, burnin=50000, sample=1000,
# #       thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.irt,"try-error")) {return()}
#     return(jags.irt)
#   }, mc.cores=2 )
# )
# save (out6, file="UNvoetenRunLastTry.RData")
load ("UNvoetenRunLastTry.RData")

chainsConv <- mcmc.list(list (out6[[1]][[1]][[1]][,-c(fix.bill+194)], out6[[2]][[1]][[1]][,-c(fix.bill+194)]))  
MixChains8 <- rbind (out6[[1]][[1]][[1]], out6[[2]][[1]][[1]])

gelman.diag (chainsConv, multivariate=FALSE) # should all be close to 1; they are
gw <- geweke.diag (MixChains8)[[1]]
mean(abs(gw)>2, na.rm=TRUE) # should be close to 0.05; it is 0.0532

# Traceplots of theta parameters
# mcmcplot (chainsConv, dir="~/Dropbox/UN/"   
#           , filename="traceplotsUNvoetenLastTry"   # <= Change name here
#           , random=10)
# Excellent mixing


Theta.voeten.last <- matrix (NA, ncol=3, nrow=ncol(MixChains8[,c(1:126)]))
for (i in 1:ncol(MixChains8[,c(1:126)])){
  Theta.voeten.last[i,] <- quantile (MixChains8[,c(1:126)][,i], prob=c(0.25,0.5,0.75))
}

# Ordered preferences according to Voeten's model, based on abstentions only as 2s
ord.voeten.last <- order (Theta.voeten.last[,2])

# To recover an order comparable to other runs
ch.ord.voeten.last <- numeric()
for (i in 1:length(y.country.order)) {
  cname <- IRT.country.order[i]
  ch.ord.voeten.last[i] <- which (y.country.order==cname)
}

y.country.order.last <- y.country.order[ch.ord.voeten.last]

# Some comparisons
plot (Theta.voeten.last[ch.ord.voeten.last,2], Theta.voeten[ch.ord,2])
plot (Theta.voeten.last[ch.ord.voeten.last,2], Theta.voeten.2[ch.ord.voeten.2,2])
plot (Theta.irt[,2], Theta.voeten.last[ch.ord.voeten.last,2])

# Same non-sensical order